
set maxvar 20000

use "data\in\estsample.dta", clear 

drop TNPE shoulder TsE NonPeakshoulder Effektpris-Effektprisshoulderelcar 
g byte PeakxEvent = Peak*EventDay

* day before 
g byte EventDay_pre = 0 
g byte EventDay_pre2 = 0 
g byte EventDay_pre3 = 0 
foreach d in 21893 21902 21937 21944 21958 21971 21979 {
  replace EventDay_pre = 1 if dato == `d' - 1 
  replace EventDay_pre2 = 1 if dato == `d' - 2 
  replace EventDay_pre3 = 1 if dato == `d' - 3 
}
forval h = 0/23 {
g byte Tpre`h'E = Treat * (hour == `h') * EventDay_pre * (1 - elcar)
label variable Tpre`h'E "`h'"
g byte Tpre`h'Eelcar = Treat * (hour == `h') * EventDay_pre * elcar
label variable Tpre`h'Eelcar "`h'"
}
* treatment day 
forval h = 0/23 {
g byte T`h'E = Treat * (hour == `h') * EventDay * (1 - elcar)
label variable T`h'E "`h'"
g byte T`h'Eelcar = Treat * (hour == `h') * EventDay * elcar
label variable T`h'Eelcar "`h'"
}
* day after 
g byte byte EventDay_post = 0 
foreach d in 21893 21902 21937 21944 21958 21971 21979 {
  replace EventDay_post = 1 if dato == `d' + 1 
}
forval h = 0/23 {
g byte Tpost`h'E = Treat * (hour == `h') * EventDay_post * (1 - elcar)
label variable Tpost`h'E "`h'"
g byte Tpost`h'Eelcar = Treat * (hour == `h') * EventDay_post * elcar
label variable Tpost`h'Eelcar "`h'"
}

forval h = 0(2)22 {
label variable Tpre`h'E "`h'"
label variable Tpre`h'Eelcar "`h'"
label variable T`h'E "`h'"
label variable T`h'Eelcar "`h'"
label variable Tpost`h'E "`h'"
label variable Tpost`h'Eelcar "`h'"
}
forval h = 1(2)23 {
label variable Tpre`h'E ""
label variable Tpre`h'Eelcar ""
label variable T`h'E ""
label variable T`h'Eelcar ""
label variable Tpost`h'E ""
label variable Tpost`h'Eelcar ""
}
/*
foreach y in TPnextdays TNPnextdays  {
g `y'elcar = `y'*elcar 
}
*/

keep lnforbruk T* *elcar hour temp* maalepktnr dato EventDay* PeakxEvent
drop TPE 
g byte weekday = dow(dato)

g byte Peak = inrange(hour,4,10)
areg lnforbruk i.dato Tpre* T0E T1E T2E T3E T4E T5E T6E T7E T8E T9E T10E T11E T12E T13E T14E T15E T16E T17E T18E T19E T20E T21E T22E T23E T?Eelcar T??Eelcar ///
  Tpost* ///
    hour hour#elcar PeakxEvent PeakxEvent#elcar /// 
  c.temp##c.temp##c.temp temp_*  c.temp#elcar ///
  /// if dato <= date("2019_12_10","YMD") + 2 ///
  , a(maalepktnr) vce(cl maalepktnr) 
est store ev2_3
estadd ysumm  

estout ev2_3 using "data\out\fig_4_results.txt", replace type ///
  keep(T*  *hour *elcar) ///
  cells(b( fmt(3)) se(par fmt(3)) _star) indicate(temp *dato*) mlabels(1 2 3 4 5 6 ) collabels(none) legend ///
  varlabels(TPE TP TNPE TNP) starlevels(* 0.10 ** 0.05 *** 0.01) ///
  stats(ymean r2 N, fmt(3 3 0 0)) //style(tex)



* Create figure from stored results
import delimited "data\out\fig_4_resultstxt", clear

drop if _n == 1 
rename (v1 v2) (variable b_se_star)
replace variable = variable[_n-1] if variable == ""
drop if strpos(b_se_star, "*") > 0 | b_se_star == "" 
keep if strpos(variable, "T") > 0 
g b = b_se_star[_n-1] if variable == variable[_n-1]
g se = b_se_star 
replace se = subinstr(se, "(", "", .)
replace se = subinstr(se, ")", "", .) 
keep if b != ""
drop b_se_star 
destring b se, replace 
g byte elcar = (strpos(variable, "elcar") > 0)
g day = "pre" if strpos(variable, "pre") > 0 
replace day = "post" if strpos(variable, "post") > 0 
replace day = "event" if day == "" 

destring variable, i(T pre elcar elcar post E Peak nextdays) g(hour)

foreach x in b {
	gen `x'_lo = `x'-1.96*se
	gen `x'_hi = `x'+1.96*se
}
*egen yhi = max(b_hi) 
*egen ylo = min(b_lo)
g yhi = 0.2 
g ylo = -0.4 
g hour1 = hour - 0.1 
g hour2 = hour 
g hour3 = hour + 0.1 

two (line b hour2 if elcar == 0, c(L) lc(black)  lp(l)) ///
  (rcap b_lo b_hi hour2 if elcar == 0, lc(black) lw(vthin)) ///
  (sc b hour3 if  elcar == 1, m(diamond)  msize(tiny) mcolor(forest_green*0.6) c(L) lc(forest_green*0.6) lp(shortdash) ) ///
  (rcap b_lo b_hi hour3 if  elcar == 1,  lc(forest_green) lw(vthin)) ///
  ///(rarea yhi ylo hour1 if inrange(hour,16,22), col(gs11%20)) /// 
  if day == "event" /// 
  , xline(16, lp(dash) lc(black)) xline(22, lp(dash) lc(black))  xlab(0(2)23, glcolor(gs13)) ylab(, glcolor(gs13)) ytitle("") xtitle("") ///xline(16) xline(21) 
  legend(order(1 "No EV" 3 "EV") col(3) pos(6)) /// 
  title("a) CPP day") /// 
  name(a, replace)

two (line b hour2 if elcar == 0, c(L) lc(black)  lp(l)) ///
  (rcap b_lo b_hi hour2 if elcar == 0, lc(black) lw(vthin)) ///
  (sc b hour3 if elcar == 1, m(diamond)  msize(tiny) mcolor(forest_green*0.6) c(L) lc(forest_green*0.6) lp(shortdash) ) ///
  (rcap b_lo b_hi hour3 if elcar == 1,  lc(forest_green) lw(vthin)) ///
  ///(rarea yhi ylo hour1 if inrange(hour,16,22), col(gs11%20)) /// 
  if day == "post" /// 
  , xline(16, lp(dash) lc(black)) xline(22, lp(dash) lc(black))  xlab(0(2)23, glcolor(gs13)) ylab(, glcolor(gs13)) ytitle("") xtitle("") ///xline(16) xline(21) 
  legend(order(1 "No EV" 3 "Large EV") col(2) pos(6)) /// 
  title("b) Day after CPP day") /// 
  name(b, replace)

grc1leg2 a b, ycommon xcommon scale(1.5) name(fig4, replace) ysize(3)
graph export "output\fig4.pdf", as(pdf) name("fig4") replace 
